%%  Assignment 5R4 for WI3097TU-p
% This is the published matlab-file

% Authors:      Martin Tichem, Bas des Tombe
% Date:         2012-06-04 
% Available on: https://code.google.com/p/wi3097tu/

% Clear Memory for Initialization
clear all;close all;clc

%%  Question 11

clear all;
fprintf('\n___________________\n### QUESTION 11 ###\n\n')

% Setting Up the Grid
xleft   = 0; xright = 10; h = 0.25;

gr.xGr  = unique([xleft:h:xright xright]);      % gridpoints
gr.dx   = diff(gr.xGr);                         % stepsizes/gridsizes
gr.xm   = 0.5*(gr.xGr(1:end-1)+gr.xGr(2:end));  % gridnodes
gr.Nx   = numel(gr.xGr);                        % number of gridpoints
gr.nx   = numel(gr.xm);                         % number of gridnodes
gr.N    = gr.xGr(2:end-1);                      % interior gridpoints, denoted n in examples, 2 boundary conditions
gr.n    = length(gr.N);                         % number of interior grid points

% Given Parameters
E       = 2e11;                                 % [N/m2] Youngs modulus
b       = 0.04;                                 % [m] width
d       = 0.2;                                  % [m] depth
I       = (1/12)*b*d^3;                         % [m4] moment of inertia
rho     = 7800;                                 % [kg/m3] mass density of the beam
g       = 9.8;                                  % [m/s2] gravitational acceleration
m       = 250;                                  % [kg] added mass
s       = 0.5;                                  % [m] width of the added mass

% Weigth
q_own   = rho*b*d*g;                            % [N/m] own weight
Q_own   = q_own*ones(length(gr.N),1);           % [N] added weight per cell
q_added = g*m*(1/s);                            % [N/m] added weight
Q_added = zeros(length(gr.N),1);                % preallocate Q_added
Q_added(gr.N>4.75 & gr.N<5.25) = q_added;       % [N/m] Amount forces applied per step
Q_added(gr.N==4.75) = (q_added/2);              % point of discontinuity
Q_added(gr.N==5.25) = (q_added/2);              % point of discontinuity
Q_tot   = Q_own + Q_added;                      % [N] own weight

% Setting Up the Matrix and Equation
A       = diag([5;6*ones(gr.n-2,1);5])+...      % matrix A
          diag(-4*ones(gr.n-1,1),-1)+...
          diag(-4*ones(gr.n-1,1),1)+...
          diag(ones(gr.n-2,1),-2)+...
          diag(ones(gr.n-2,1),2); 
w       = [0;((E*I)/(h^4))*sparse(A)\Q_tot;0];  % define w for entire grid
x_max   = gr.xGr(w==max(w));                    % find the position of the max deflection
fprintf('The maximum deflection of the beam is:\t %1.7e [m]\nMax deflection is at x is: \t\t %1.1d [m]\n',max(w),x_max)

% Table
fprintf('\n---TABLE---\nx_pos \t\t Deflection \n')
for i = 1:4:length(gr.N)+4
    fprintf('%1.1e \t %1.7e [m]\n',(i-1)/4,w(i))% position, deflection
end

%%  Question 12

clear all;
fprintf('\n___________________\n### QUESTION 12 ###\n\n')

% Setting Up the Grid
xleft   = 0; xright = 10; h = 0.05;

gr.xGr  = unique([xleft:h:xright xright]);      % gridpoints
gr.dx   = diff(gr.xGr);                         % stepsizes/gridsizes
gr.xm   = 0.5*(gr.xGr(1:end-1)+gr.xGr(2:end));  % gridnodes
gr.Nx   = numel(gr.xGr);                        % number of gridpoints
gr.nx   = numel(gr.xm);                         % number of gridnodes
gr.N    = gr.xGr(2:end-1);                      % interior gridpoints, denoted n in examples, 2 boundary conditions
gr.n    = length(gr.N);                         % number of interior grid points

% Given Parameters
E       = 2e11;                                 % [N/m2] Youngs modulus
b       = 0.04;                                 % [m] width
d       = 0.2;                                  % [m] depth
I       = (1/12)*b*d^3;                         % [m4] moment of inertia
rho     = 7800;                                 % [kg/m3] mass density of the beam
g       = 9.8;                                  % [m/s2] gravitational acceleration
m       = 250;                                  % [kg] added mass
s       = 0.5;                                  % [m] width of the added mass

% Weigth
q_own   = rho*b*d*g;                            % [N/m] own weight
Q_own   = q_own*ones(length(gr.N),1);           % [N] added weight per cell
q_added = g*m*(1/s);                            % [N/m] added weight
Q_added = zeros(length(gr.N),1);                % preallocate Q_added
Q_added(gr.N>4.75 & gr.N<5.25) = q_added;       % [N/m] Amount forces applied per step
Q_added(gr.N==4.75) = (q_added/2);              % point of discontinuity
Q_added(gr.N==5.25) = (q_added/2);              % point of discontinuity
Q_tot   = Q_own + Q_added;                      % [N] total weight

% Setting Up the Matrix and Equation
A       = diag([5;6*ones(gr.n-2,1);5])+...      % matrix A
          diag(-4*ones(gr.n-1,1),-1)+...
          diag(-4*ones(gr.n-1,1),1)+...
          diag(ones(gr.n-2,1),-2)+...
          diag(ones(gr.n-2,1),2); 
w1      = [0;((E*I)/(h^4))*sparse(A)\Q_own;0];  % define w1 for entire grid
w2      = [0;((E*I)/(h^4))*sparse(A)\Q_added;0];% define w2 for entire grid

% Plot Function
g       = figure();
plot(gr.xGr,w1,'+',gr.xGr,w2,'o');
title('Question 12 - Beam Deflection');
xlabel('Position [m]');
ylabel('Deflection [m]');
set(gca,'ydir','reverse');
legend('Own Weight','Added Weight','Location','Southeast');

% Return the values at midpoint
fprintf('Deflection due to own   weight at x=5:\t %1.7e [m]\nDeflection due to added weight at x=5:\t %1.7e [m]\n',max(w1),max(w2))

%% Question 13

clear all;
fprintf('\n___________________\n### QUESTION 13 ###\n\n')

% Setting Up the Grid
xleft   = 0; xright = 10; h = 0.05;

gr.xGr  = unique([xleft:h:xright xright]);      % gridpoints
gr.dx   = diff(gr.xGr);                         % stepsizes/gridsizes
gr.xm   = 0.5*(gr.xGr(1:end-1)+gr.xGr(2:end));  % gridnodes
gr.Nx   = numel(gr.xGr);                        % number of gridpoints
gr.nx   = numel(gr.xm);                         % number of gridnodes
gr.N    = gr.xGr(2:end-1);                      % interior gridpoints, denoted n in examples, 2 boundary conditions
gr.n    = length(gr.N);                         % number of interior grid points

% Given Parameters
E       = 2e11;                                 % [N/m2] Youngs modulus
b       = 0.04;                                 % [m] width
d       = 0.2;                                  % [m] depth
I       = (1/12)*b*d^3;                         % [m4] moment of inertia
rho     = 7800;                                 % [kg/m3] mass density of the beam
g       = 9.8;                                  % [m/s2] gravitational acceleration
m       = 250;                                  % [kg] added mass
s       = 0.5;                                  % [m] width of the added mass

% Weigth
q_own   = rho*b*d*g;                            % [N/m] own weight
Q_own   = q_own*ones(length(gr.N),1);           % [N] added weight per cell
q_added = g*m*(1/s);                            % [N/m] added weight
Q_added = zeros(length(gr.N),1);                % preallocate Q_added
Q_added(gr.N>4.75 & gr.N<5.25) = q_added;       % [N/m] Amount forces applied per step
Q_added(gr.N==4.75) = (q_added/2);              % point of discontinuity
Q_added(gr.N==5.25) = (q_added/2);              % point of discontinuity

% Find Deflection Due to Own Weight
A       = diag([5;6*ones(gr.n-2,1);5])+...      % matrix A
          diag(-4*ones(gr.n-1,1),-1)+...
          diag(-4*ones(gr.n-1,1),1)+...
          diag(ones(gr.n-2,1),-2)+...
          diag(ones(gr.n-2,1),2); 
w       = [0;((E*I)/(h^4))*sparse(A)\Q_own;0];  % define w for entire grid

% Increase the Mass With Increment 'inc' Until w_max == 2*w(q_own)

w_ref   = max(w);    % reference value for deflection
w_max   = 0;         % initiate max deflection                          
m       = 350;       % initial mass
inc     = 1;         % increment for mass

while w_max < 2*w_ref 
    m = m + inc;                                    % increment
    
    if m > 1000                                     % break in case routine fails
        break;
    end
    
    q_added = g*m*(1/s);                            % [N/m] added weight
    Q_added = zeros(length(gr.N),1);                % preallocate Q_added
    Q_added(gr.N>4.75 & gr.N<5.25) = q_added;       % [N] Amount forces applied per step
    Q_added(gr.N==4.75) = (q_added/2);              % point of discontinuity
    Q_added(gr.N==5.25) = (q_added/2);              % point of discontinuity
        
    Q_tot   = Q_own + Q_added;                      % [N] own weight per cell        
    w       = [0;((E*I)/(h^4))*sparse(A)\Q_tot;0];  % define w for entire grid
    w_max   = max(w);
end 
fprintf('The maximum allowable added mass is: \t %1.1d [kg]\n',m)

%% Question 14

clear all;
fprintf('\n___________________\n### QUESTION 14 ###\n\n')

% Setting Up the Grid
xleft   = 0; xright = 10; h = .05;

gr.xGr  = unique([xleft:h:xright xright]);      % gridpoints
gr.dx   = diff(gr.xGr);                         % stepsizes/gridsizes
gr.xm   = 0.5*(gr.xGr(1:end-1)+gr.xGr(2:end));  % gridnodes
gr.Nx   = numel(gr.xGr);                        % number of gridpoints
gr.nx   = numel(gr.xm);                         % number of gridnodes
gr.N    = gr.xGr(2:end-1);                      % interior gridpoints, denoted n in examples, 2 boundary conditions
gr.n    = length(gr.N);                         % number of interior grid points

% Given Parameters
E       = 2e11;                                 % [N/m2] Youngs modulus
b       = 0.04;                                 % [m] width
d       = 0.2;                                  % [m] depth
I       = (1/12)*b*d^3;                         % [m4] moment of inertia
rho     = 7800;                                 % [kg/m3] mass density of the beam
g       = 9.8;                                  % [m/s2] gravitational acceleration
m       = 250;                                  % [kg] added mass
s       = 0.5;                                  % [m] width of the added mass

% Weigth
q_own   = rho*b*d*g;                            % [N/m] own weight
Q_own   = q_own*ones(length(gr.N),1);           % [N] added weight per cell
q_added = g*m*(1/s);                            % [N/m] added weight
Q_added = zeros(length(gr.N),1);                % preallocate Q_added
Q_added(gr.N>4.75 & gr.N<5.25) = q_added;       % [N/m] Amount forces applied per step
Q_added(gr.N==4.75) = (q_added/2);
Q_added(gr.N==5.25) = (q_added/2);

% Calculating Value of Am by Intergrating Eq(10)
syms x                                          % symbolic x
f   = sin(pi*(x-4.75)/s);                       % weight distribution function (10)
Am  = double(m*g/(int(f,4.75,5.25)*g));         % solving for Am

q_added2 = Am*g*sin(pi*(gr.N-4.75)/s);          % [N/m] discrete distribution of q_m                  
Q_added2 = zeros(length(gr.N),1);               % preallocate Q_added2
Q_added2(gr.N>4.75 & gr.N<5.25) ...             % [N/m] distrubution vector of q_m 
         = q_added2(gr.N>4.75 & gr.N<5.25);   
  
% Q_added2(gr.N==4.75) = (q_added2(gr.N==4.75)/2);
% Q_added2(gr.N==5.25) = (q_added2(gr.N==5.25)/2);

% Drawing the load
g       = figure();
plot(gr.N,Q_added,gr.N,Q_added2)
title('Question 14 - Load Distribution');
xlabel('Position [m]');
ylabel('Load [N/m]');
legend('Even Distributed Load','Curved Load','Location','Northeast');
set(gca,'YTick',-500:1000:6000)
axis([4.5 5.5 -2000 9000])

fprintf('The calculated value of Am is: \t\t %1.7d []\n',Am)
%% Question 15

clear all;
fprintf('\n___________________\n### QUESTION 15 ###\n\n')

% Setting Up the Grid
xleft   = 0; xright = 10; h = .05;

gr.xGr  = unique([xleft:h:xright xright]);      % gridpoints
gr.dx   = diff(gr.xGr);                         % stepsizes/gridsizes
gr.xm   = 0.5*(gr.xGr(1:end-1)+gr.xGr(2:end));  % gridnodes
gr.Nx   = numel(gr.xGr);                        % number of gridpoints
gr.nx   = numel(gr.xm);                         % number of gridnodes
gr.N    = gr.xGr(2:end-1);                      % interior gridpoints, denoted n in examples, 2 boundary conditions
gr.n    = length(gr.N);                         % number of interior grid points

% Given Parameters
E       = 2e11;                                 % [N/m2] Youngs modulus
b       = 0.04;                                 % [m] width
d       = 0.2;                                  % [m] depth
I       = (1/12)*b*d^3;                         % [m4] moment of inertia
rho     = 7800;                                 % [kg/m3] mass density of the beam
g       = 9.8;                                  % [m/s2] gravitational acceleration
m       = 250;                                  % [kg] added mass
s       = 0.5;                                  % [m] width of the added mass

% Weigth
q_own   = rho*b*d*g;                            % [N/m] own weight
Q_own   = q_own*ones(length(gr.N),1);           % [N] added weight per cell
q_added = g*m*(1/s);                            % [N/m] added weight
Q_added = zeros(length(gr.N),1);                % preallocate Q_added
Q_added(gr.N>4.75 & gr.N<5.25) = q_added;       % [N/m] Amount forces applied per step
Q_added(gr.N==4.75) = (q_added/2);              % point of discontinuity
Q_added(gr.N==5.25) = (q_added/2);              % point of discontinuity

% Calculating Value of Am by Intergrating Eq(10)
syms x                                          % symbolic x
f   = sin(pi*(x-4.75)/s);                       % weight distribution function (10)
Am  = double(m*g/(int(f,4.75,5.25)*g));         % solving for Am

q_added2 = Am*g*sin(pi*(gr.N-4.75)/s);          % [N/m] discrete distribution of q_m                  
Q_added2 = zeros(length(gr.N),1);               % preallocate Q_added2
Q_added2(gr.N>4.75 & gr.N<5.25) ...             % [N/m] distrubution vector of q_m 
         = q_added2(gr.N>4.75 & gr.N<5.25);   

% Setting Up the Matrix and Equation
A       = diag([5;6*ones(gr.n-2,1);5])+...       % matrix A
          diag(-4*ones(gr.n-1,1),-1)+...
          diag(-4*ones(gr.n-1,1),1)+...
          diag(ones(gr.n-2,1),-2)+...
          diag(ones(gr.n-2,1),2); 
w1      = [0;((E*I)/(h^4))*sparse(A)\Q_own;0];   % define w1 for entire grid
w2      = [0;((E*I)/(h^4))*sparse(A)\Q_added;0]; % define w2 for entire grid
w3      = [0;((E*I)/(h^4))*sparse(A)\Q_added2;0];% define w3 for entire grid

% Plot Function
g       = figure();
plot(gr.xGr,w2,'+',gr.xGr,w3,'o');
title('Question 15 - Beam Deflection');
xlabel('Position [m]');
ylabel('Deflection [m]');
set(gca,'ydir','reverse');
legend('Even Distributed Load','Curved Load','Location','North');

% Return the values at midpoint
fprintf('Max deflection even distributed load:\t %1.7e [m]\nMax deflection curved load: \t\t %1.7e [m]\n\n\n',max(w2),max(w3))
